# 26_utas2017.R
# Purpose: Ideological Extremism and Political Participation in Japan
# Created: 2021-3-14 Taka-aki Asano
# Last Modified: 2021-10-10

# package
require("ltm")
require("plink")
require("readr")
require("dplyr")
require("tidyr")
require("psych")


# Voters
## read dataset
voter2017 <- read_csv(
  "http://www.masaki.j.u-tokyo.ac.jp/utas/2017UTASV20200326.csv", 
  locale = locale(encoding = "CP932")
)

## specify variables
variables2017 <- c("Q22_1", "Q23_1", "Q23_2", "Q23_3", "Q23_4", 
                   "Q23_5", "Q23_6", "Q23_7", "Q23_8", "Q23_9", 
                   "Q23_10", "Q23_11", "Q23_12", "Q23_13", "Q23_14", 
                   "Q23_15", "Q23_16", "Q23_17", "Q24_1", "Q24_2", 
                   "Q24_3", "Q24_4", "Q24_5", "Q25_1", "Q25_2", 
                   "Q25_3", "Q25_4", "Q25_5", "Q26_1", "Q27")

## handle missing data
for(i in variables2017) {
  n <- voter2017[,i] == 66 | voter2017[,i] == 99
  voter2017[n, i] <- NA
}
voter2017 <- voter2017[rowSums(voter2017[,variables2017], na.rm = TRUE) != 0,]

## estimate item parameters
irt2017_voter <- ltm::grm(voter2017[,variables2017], IRT.param = FALSE)
irt2017_voter ## viewing

## estimate voters ideology
score2017 <- ltm::factor.scores(irt2017_voter, voter2017[,variables2017])


# rescale
## common item
common2009 <- c("Q013801", "Q013802", "Q013803", "Q013804", "Q013806", 
                "Q013809", "Q013811", "Q013812", "Q013817", "Q013818")
common2017 <- c("Q26_1", "Q23_1", "Q23_4", "Q23_2",  "Q23_3", 
                "Q23_6", "Q23_7", "Q23_8", "Q23_12", "Q23_11")
index <- data.frame(
  group1 = match(common2009, colnames(voter2009[,variables2009])), 
  group2 = match(common2017, colnames(voter2017[,variables2017]))
)
index

## list of parameters
### 2009
res2009 <- coef(irt2009_voter)
res2009 <- res2009[,c(5,1:4)]
colnames(res2009) <- c("a","b1","b2","b3","b4")
### 2017
res2017 <- coef(irt2017_voter)
res2017 <- res2017[,c(5,1:4)]
colnames(res2017) <- c("a","b1","b2","b3","b4")
### merge
pm <- list(res2009, res2017)

## list of scale
rescat <- list(rep(5, nrow(res2009)), 
               rep(5, nrow(res2017)))

## as.irt.pars
pm2009 <- as.poly.mod(n = nrow(res2009), model = "grm", 
                      items = 1:nrow(res2009))
pm2017 <- as.poly.mod(n = nrow(res2017), model = "grm", 
                      items = 1:nrow(res2017))
res <- as.irt.pars(pm, common = index, cat = rescat, 
                   poly.mod = list(pm2009, pm2017), 
                   location = FALSE)

## list of voters' position
ideology <- list(score2009$score.dat$z1, score2017$score.dat$z1)

## rescale
plink_out <- plink(res, common = index, base.grp = 1, 
                   rescale = "MS", ability = ideology)
summary(plink_out)

## estimate voters ideology (rescale)
score2017_plink <- link.ability(plink_out)$group2
voter2017$Ideology <- -1 * score2017_plink

## Median
median(voter2017$Ideology)
by(voter2017$Ideology, voter2017$Q2, describe)
